Selection for stress tolerance and longevity in Drosophila melanogaster have strong impacts on microbiome profiles

There is experimental evidence that microbiomes have a strong influence on a range of host traits. Understanding the basis and importance of symbiosis between host and associated microorganisms is a rapidly developing research field, and we still lack a mechanistic understanding of ecological and genetic pressures affecting host-microbiome associations. Here Drosophila melanogaster lines from a large-scale artificial selection experiment were used to investigate whether the microbiota differ in lines selected for different stress resistance traits and longevity. Following multiple generations of artificial selection all selection regimes and corresponding controls had their microbiomes assessed. The microbiome was interrogated based on 16S rRNA sequencing. We found that the microbiome of flies from the different selection regimes differed markedly from that of the unselected control regime, and microbial diversity was consistently higher in selected relative to control regimes. Several common Drosophila bacterial species showed differentially abundance in the different selection regimes despite flies being exposed to similar environmental conditions for two generations prior to assessment. Our findings provide strong evidence for symbiosis between host and microbiomes but we cannot reveal whether the interactions are adaptive, nor whether they are caused by genetic or ecological factors.

of environmental-stress or life-history traits selected for increased cold-shock resistance, desiccation resistance, heatshock resistance, all measured as survival proportion, starvation resistance, heat-knockdown resistance and longevity respectively measured as survival duration in days, survival after exposure to heat stress in minutes, and LT50 in hours.
Phenotypic traits were recorded for flies from five replicate lines for each of the selection regimes (CS: cold resistance; DS: desiccation resistance; HS: heat shock resistance; SS: starvation resistance; KS: heat-knock down resistance; LS: longevity) and for flies from five unselected control lines (UC) and mean phenotypic values were estimated per replication line.Observed, Shannon's, and Faith's Phylogenetic alpha-diversity metrics were estimated per replication line using taxonomy filtered, pruned and rarefied amplicon sequence variant (ASV) data.For each selection regime, alpha diversity of the five UC lines and the five replication lines under selection were plotted against their mean phenotypic trait values and the Pearson´s correlation coefficient (r) between alpha diversity and phenotypic trait values were estimated.Colour code applied indicates affiliation with the selection regime.Dashed black lines represents the "line of best fit", minimizing vertical distances between the data points and the line itself, representing the linear trends between alpha diversity and selection responsiveness of phenotypic traits.Selection for increased stress resistance or longevity was performed every other generation and populations were selected for 11-27 generations (depending on the selection regime) using a common garden design.Prior to phenotypic assessment, flies from all regimes were reared at benign conditions (like those at which UC flies were reared across all generations) for two generations (for details see Bubliy and Loeschcke (2005) and description in Material and Methods).Flies used to characterize gut microbial alpha diversity were collected from the same selection experiment but after 14-32 generations of selection depending on the selection regime (i.e., 32 selection events for HS, CS, KS, DS and UC; 37 selection events for SS; 14 events of selection for LS; see Materials and Methods for more details), following benign conditions for two generations prior to sample collection for gut microbiome characterization.Unfortunately, phenotypic assessment of selection responses was not performed in the later generations where flies were harvested for microbiome analysis, hence, correlations between microbial diversity and phenotypic responsiveness to selection could not be investigated within the same individual and should only be considered as proxy.baseMean: Mean abundance of the normalized count values over all samples included in the comparison; log2FC: log2-fold change represents the effect size estimate reported on a logarithmic scale to base 2, indicating how much the ASV abundance has changed due to selection for heat knockdown resistance in comparison to the control lines; P adj: adjusted P-values using the Benjamini-Hochberg procedure.
Table S8.Differential abundance analysis comparing lines selected for longevity (LS) with unselected control lines (UC) at genus level.Taxonomies were agglomerated to genus level prior to analysis.ASVs are ordered according to adjusted P-values (P adj) from lowest to highest.ASVs were considered differentially abundant for adjusted Pvalues ≤ 0.01 and an FDR cut-off of 5%.Significance threshold (P adj = 0.01) is marked by a dashed line.

Figure S2 .
Figure S2.Scatter plots showing the relationship between gut microbial alpha diversity (panel A: observed alpha diversity; panel B: Shannon´s diversity index; panel C: Faith´s phylogenetic diversity index) and phenotypic measures

Table S1 .
Non-parametric pairwise test statistics (Wilcoxon Rank Sum test with Benjamini-Hochberg multiple hypothesis test correction) for observed -diversity (species richness measured as number of observed ASVs).Rows in bold indicate statistically significant differences between pairs of selection regimes.
∆ difference -diversity between pairs of selection regimes.

Table S2 .
Non-parametric pairwise test statistics (Wilcoxon Rank Sum test with Benjamini-Hochberg multiple hypothesis test correction) for Shannon's -diversity index.Rows in bold indicate statistically significant differences between pairs of selection regimes.
∆ difference -diversity between pairs of selection regimes.

Table S3 .
Parametric pairwise test statistics (Wilcoxon Rank Sum test with Benjamini-Hochberg multiple hypothesis test correction) for Faith´s phylogenetic diversity.Rows in bold indicate statistically significant differences between pairs of selection regimes.

Table S4 .
Differential abundance analysis comparing lines selected for cold-shock resistance (CS) to the unselected control lines (UC) at the genus level.Taxonomies were agglomerated to genus level prior to analysis.ASVs are ordered according to adjusted P-values (P adj) from lowest to highest.ASVs were considered differentially abundant for adjusted P-values ≤ 0.01 and an FDR cut-off of 5%.Significance threshold (P adj = 0.01) is marked by a dashed line.
baseMean: Mean abundance of the normalized count values over all samples included in the comparison; log2FC: log2-fold change represents the effect size estimate reported on a logarithmic scale to base 2, indicating how much the ASV abundance has changed due to selection for cold-shock resistance in comparison to the control lines; P adj: adjusted P-values using the Benjamini-Hochberg procedure.

Table S5 .
Differential abundance analysis comparing lines selected for desiccation resistance (DS) to the unselected control lines (UC) at the genus level.Taxonomies were agglomerated to genus level prior to analysis.ASVs are ordered according to adjusted P-values (P adj) from lowest to highest.ASVs were considered differentially abundant for adjusted P-values ≤ 0.01 and an FDR cut-off of 5%.Significance threshold (P adj = 0.01) is marked by a dashed line.
baseMean: Mean abundance of the normalized count values over all samples included in the comparison; log2FC: log2-fold change represents the effect size estimate reported on a logarithmic scale to base 2, indicating how much the ASV abundance has changed due to selection for desiccation resistance in comparison to the control lines; P adj: adjusted P-values using the Benjamini-Hochberg procedure.

Table S6 .
Differential abundance analysis comparing lines selected for heat-shock resistance (HS) to the unselected control lines (UC) at the genus level.Taxonomies were agglomerated to genus level prior to analysis.ASVs are ordered according to adjusted P-values (P adj) from lowest to highest.ASVs were considered differentially abundant for adjusted P-values ≤ 0.01 and an FDR cut-off of 5%.Significance threshold (P adj = 0.01) is marked by a dashed line.
baseMean: Mean abundance of the normalized count values over all samples included in the comparison; log2FC: log2-fold change represents the effect size estimate reported on a logarithmic scale to base 2, indicating how much the ASV abundance has changed due to selection for heat-shock resistance in comparison to the control lines; P adj: adjusted P-values using the Benjamini-Hochberg procedure.

Table S7 .
Differential abundance analysis comparing lines selected for heat knockdown resistance (KS) to the unselected control lines (UC) at the genus level.Taxonomies were agglomerated to genus level prior to analysis.ASVs are ordered according to adjusted P-values (P adj) from lowest to highest.ASVs were considered differentially abundant for adjusted P-values ≤ 0.01 and an FDR cut-off of 5%.Significance threshold (P adj = 0.01) is marked by a dashed line.

Table S9 .
Differential abundance analysis comparing lines selected for starvation resistance (SS) to the unselected control lines (UC) at the genus level.Taxonomies were agglomerated to genus level prior to analysis.ASVs are ordered according to adjusted P-values (P adj) from lowest to highest.ASVs were considered differentially abundant for adjusted P-values ≤ 0.01 and an FDR cut-off of 5%.Significance threshold (P adj = 0.01) is marked by a dashed line.Mean abundance of the normalized count values over all samples included in the comparison; log2FC: log2-fold change represents the effect size estimate reported on a logarithmic scale to base 2, indicating how much the ASV abundance has changed due to selection for starvation resistance in comparison to the control lines; P adj: adjusted P-values using the Benjamini-Hochberg procedure.